###############################################################
# REPLICATE MAIN RESULTS REMOVING ONE REGION AT A TIME: COUNT #
###############################################################
# Author: Kasia Nalewajko
# First created: 3 December 2023
# Replicated: 17 June 2024

rm(list = ls())

# LOAD PACKAGES -----------------------------------------------------------

if (!require("dplyr")) install.packages("dplyr")
if (!require("fixest")) install.packages("fixest")
if (!require("ggplot2")) install.packages("ggplot2")

# LOAD DATA ---------------------------------------------------------------

load("./00 SUBMITTED/00 APSR final/04 replication_files/01 data/main.Rda")
  
# RUN MODEL IN A LOOP -----------------------------------------------------

resregions <- unique(main$resregion2)
resregions <- resregions[resregions != "divided"]
resregions <- sort(resregions)

# Create empty containers
removed_resregion <- rep(NA, 1*length(resregions))
estimate <- rep(NA, 1*length(resregions))
std_error <- rep(NA, 1*length(resregions))
models_list <- list()

# Loop over resregions
for(i in 1:length(resregions)){
  
  models_list[[i]] <- fixest::feols(log(all_jewish_victims+1) ~ log(pop1936+1) + synagogues_dummy + catholic_churches + log(collabos_antijew1000+1) + log(DHI_milipol_sum1942+1) + log(FRANCISTE_36_1_pct+1) + log(ActionFrancaise_19_pct+1) + turnout_36_1_pct + log(nuance_extremedroite_36_1_pct+1) + log(nuance_centredroit_36_1_pct+1) + log(nuance_droite_36_1_pct+1) + log(nuance_centregauche_36_1_pct+1) + log(nuance_gauche_36_1_pct+1) + log(nuance_extremegauche_36_1_pct+1) + area_sqkm  + longitude + longitudesq + latitude + latitudesq | zone5 + as.factor(ar_name)
                                    | log((FFIFFCgentile*1000/pop1936)+1) ~ log(mpf1000+1),
                                    cluster = ~c(id_bureau, resregion2), 
                                    data = subset(main, resregion2 != resregions[i]))
  
  removed_resregion[[i]] <- resregions[i]
  estimate[[i]] <- models_list[[i]][["coefficients"]][["fit_log((FFIFFCgentile * 1000/pop1936) + 1)"]]
  std_error[[i]] <- models_list[[i]][["se"]][["fit_log((FFIFFCgentile * 1000/pop1936) + 1)"]]
  
}

models_df <- data.frame(
  estimate = estimate,
  std_error = std_error,
  removed_resregion = removed_resregion,
  ci_hi_95 = estimate + 1.96 * std_error,
  ci_lo_95 = estimate - 1.96 * std_error,
  ci_hi_90 = estimate + 1.645 * std_error,
  ci_lo_90 = estimate - 1.645 * std_error)

models_df$removed_resregion <- paste0("Without ", models_df$removed_resregion)

# PLOT -----------------------------------------------------

models_df %>%
  ggplot(aes(
    x = removed_resregion,
    y = estimate
  )) +
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") +
  geom_linerange(aes(ymin = ci_lo_90, ymax = ci_hi_90), 
                 position = position_dodge(width = 0.1),
                 size = 0.5) +
  geom_point(size = 1.75
             , position = position_dodge(width = 0.1)
  ) +
  theme_bw() +
  labs(
    x = "Removed region",
    y = "Estimates"
  ) +
  coord_flip()

# EXPORT -----------------------------------------------------

ggsave(
  "C1_removing_regions_count.png",
  plot = last_plot(),
  path = "./00 SUBMITTED/00 APSR final/02 supplementary_materials/figures",
  width = 20,
  height = 24,
  units = "cm",
  dpi = 300
)



